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I. 



INTRODUCTION 



Considerable interest has recently been shown in trap- 
ped waves travelling along the boundaries of continents. A 
"waveguide" effect exists over the continental shelf. That 
is, wave energy is confined (essentially by refraction) to 
the continental shelf. Two general types of these waves 
exist : 

A. Edgewaves which are characterized by wavelengths of hun- 
dreds of kilometers (km) and periods of hours (almost always 
less than a pendulum day) . 

B. Shelf waves, which are generally even longer, have pe- 
riods greater than one pendulum day, and travel southward 
along the west coast of an ocean in the northern hemisphere 
(as do Kelvin waves) . 

STOKES (1846) showed that such a wave guide effect ex- 
ists over a uniformly sloping beach or continental shelf, 
with the amplitude of the gravity waves decaying exponen- 
tially to seaward. URSELL (1952) showed that Stoke' s edge- 
waves were the fundamental mode of a family of waves order- 
ed by the number of modes parallel to the coast. REID(1958) 
studied long waves on uniformly sloping shelves of infinite 
width, including the effect of the Coriolis force. Reid 
showed that the sea surface may react as an "inverse baro- 
meter" and that atmospheric pressure systems may be a driv- 
ing force for edgewaves. He found that the Coriolis force 
could cause the wave period to vary from 46% less than to 
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86% greater than that for the non-rotating case, depending 
on the direction of travel. A new quasigeos trophic wave is 
now possible, analogous to a Kelvin wave, having.no small 
scale counterpart. 

ROBINSON (1964) ini tiated a study of the continental 
shelf wave and studied the data of HAMON (1962, 1963) re- 
lating tidal and barometric conditions at several stations 
on the eastern and western coasts of Australia. In this 
model the continental shelf ends abruptly, at which point 
the depth becomes infinite. He found that an inverse baro- 
meter effect was exhibited but that the propagated shelf 
waves had a celerity double that of his calculations for 
the western boundary. MOOERS and SMITH (1968) studied the 
relation of sea level and barometric conditions along the 
Oregon coast for a period of nearly one year. Their statis- 
tical results show a barometric factor of -1.2 cm/mb and 
predominant sea level oscillations of 0.1 and 0.35 cycles 
per day in the summer. They conclude that a shelf wave of 
period three days is travelling north. MYSAK and HAMON 
(1969) found shelf waves off the coast of North Carolina, 
but found no coupling between the sea surface and air pres- 
sure in the frequency range 0-0.5 cpd. ADAMS and BUCHWALD 
(1969) show that an equally suitable driving force for 
shelf waves is the longshore component of the geos trophic 
wind. This may account for the exaggerated frequency re- 
sponse of the sea level on the east coast observed by 
Hamon. 
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MYSAK (1967, 1968) extended Robinson ' s work, and discus- 
sed the effect of a continental shelf of finite width on the 
frequency of Hamon 1 s Australian waves. His theoretical so- 
lutions correspond more closely to the observations, al- 
though he still cannot account for the extremely low read- 
ings along the eastern boundary. He attributes the discrep- 
ancy mainly to the presence of stratified water and currents 
in the deep water beyond the continental shelf. A signifi- 
cant discrepancy exists between the dispersion relation and 
that for waves over an infinitely wide continental shelf. 

This paper is a study of the effects of a continental 
slope and finite ocean depth upon the present one-slope 
models of MYSAK (1968) . A sharp discontinuity in the depth 
of water beyond the continental shelf is not a common occur- 
rence in the world ocean. It is interesting to study the 
two-slope situation where a gently sloping continental shelf 
(slope, s < 0.002) and steeper continental slope (s 0.05) 
form a transition zone between the coast-line and deep 
water. Three parameters: the slope of the continental 

shelf, the slope of the continental slope, and the depth of 
the deep water should have possible effects on shelf waves. 
These are investigated below. 
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II. ANALYSIS 



h^= s^x(x ^ W^) 




irrrr 



V S 2 X + W l (s l- S 2> 



h^= constant 



(W x £ x £ W 2 + Vi^) 



I 



II 



III 




Fig. 1. Profile of the two-slope model. 

The model is characterized by a gradually sloping con- 
tinental shelf of finite width (Region I) adjoining a steep- 
er continental slope (Region II) which terminates in water 
of uniform depth (Region III) . A representative slope and 
width of the continental shelf of .002 and 100 km (Mysak, 
1968a) are used below. A representative depth of the 

deep ocean is 5000m and is the greatest depth of Region III. 
Three slopes will be used for the continental slope: .03, 

.05 and .08, with .05 used as the standard for comparison 
with Mysak" s model. 

The shallow water equations are used: 




d/dx(hu) + 5/dy(hv) 4 - 6£/dt = 0 



( 2 ) 
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where (u, v) are the (x, y) velocity components and C 
is the free surface height. 

Consider a wave, moving in the y direction, speci- 
fied by 




i i k ( X ) e i < ot - m y> 

V k (x)e i(<Tt - m y> 

n k (x)e i(fft - m ^ ) 



(3) 



where k denotes the region. 
Using (1) and (3) , u and v are 



u. 



v. 



-^~2 

f -a 



- y - - V (crmr? 
f -<J 



9J7 x i(gt-my) 
dx } k e 



d r? i(gt-my) 
5y ; k e 



(4) 

(5) 



Eliminating u, v from (2) in Region I gives the equa- 



tion 



h l 7? l + S 1 T7 1 + ( 



_2 ^2 s,fm 0 

i- - V ) ”i ■ 0 



( 6 ) 



After making the substitutions 



h i - 
p i = 



s^x 

2 _2 
g -f 

gs. 



fm 

g 



(6) can be written 



xT^" + r^' + (p 1 -m 2 x) T? 1 = 0 



(7) 
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This equation in has as its solution (REID, 1958) 

- z./2 

r?i = CA 1 Fa 1 (z 1 ) + B 1 Ga 1 (z 1 )}e (8) 



where Fa^(z^) is Rummer's function, and Ga^(z^) is a sec- 
ond (independent) solution (SLATER, 1960) : 

^ \ . . a(a+l)z 2 , a(a+l) (a+2)z 3 , 

Fa(z) =l + az+ — ^ J + — » L + • • • 

+ a (a+1) ♦ ♦ • (a+n-1) z n + 

(ni ) 2 



and 
Ga (z) 

+ 



= Fa (z) In z + az (— - 2)+ a (a+1) z 2 (^a+lf - ^ * * * 

a (a+1) • • • (a+n-1) z n (— H — + • • • + —• - .. - 2-1- •••—)+ ••• 
v ' v ' 'a a+1 a+n-1 n 



( 10 ) 



Here, a^ = 
of integration, 
approaches zero 



V 



1 



A l Fa l 



m-Pf 

’ 2m ' 



z^ = 2mx, and A^ , are constants 



Since Ga^(z^) approaches infinity as z 



, B, must equal zero: 



(z i )e 



~ z j/ 2 



( 11 ) 



Substituting (11) into (4) gives 



z ]_/2 dFa (z ) 

U 1 = ' A d( f+CT )Pa 1 ( z 1 )-2< r -a57— }=° < 12 > 



a -f 



Similarly, in Region II, 

— z /2 

V 2 = e {A 2 Fa 2 (z 2 ) + B 2 Ga 2 (z 2 )} (13) 
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and 



z 2 /2 

TT _ _ iqme 

2 2 2 
or -f 



{A 2 [(a-f)Fa 2 (z 2 )- 2a dz 



dGa 2 (z 2 ) - 



+ B 2 [(a-f)Ga 2 (z 2 )- 2a ~ 



dGa 2 (z 2 ^ ^ 



(14) 



where 



z 2 = 2m? 2 , a 2 



m-p, 

2m 



? 2 = and ?2 



g 2 -f 2 

gs„ 



fm 

a 



In Region III the counterpart of (6) is 



h 3”3" + ( £_ 5^ _ - h 3 m2 )n 3 = 0 



(15) 



Since r\ must be bounded for large x. 



•a 

ri ^ = A^e where 



r o 2 -2 -a 1 

= L m - -35T ] 2 and 



(16) 



igA 



U 3 2 2 

J cr -f 



3 -£(x-W -,-W 2 ) 



(fm+cy^) e 



(17) 



There are now equations defining T) and U in each re- 
gion. The next step is to patch together the solutions for 
f) and U at the points x = W^, and x = (W^+W 2 ) thus eliminat- 
ing the constants A^ 



, . The patching conditions are 
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za[ 


- C « ] 2=° 


(surface height continuity) 


(18) 


Chu3 1 


3 

= [hu] 2 = 0 


(normal flux continuity) 


(19) 


where 

[ ] 3 
k 


V [ 


■*k * 





The following abbreviations will be used; 



F. = Fa.(z.) G. = Ga.(z.) 
1 3 3 3 3 3 



The subscript, 1, refers to the solution for Region I 
(continental shelf) where it joins Region II (continental 
slope). The subscript, 2, refers to the Region II where it 
joins Region I. The subscript, 3, refers to the continental 
slope where it joins Region III, the flat bottom. The func- 
tions subscripted 3 have the same form as those subscripted 
2 with the exception of the variable, z^» which is deter- 
mined by the distance from the origin. 

Using (18) and (19) between Regions I and II and set- 
ting A^ = 1, the constants of integration A^ and can be 
solved for? 



A 2 = 



G 2 F l' F 1 G 2 
G 2 F 2 " F 2 G 2 



B 2 = 



F 1 F 2 ' ~ F 2 F l' 
G 2 F 2 ' ~ F 2 G 2 ' 



( 20 ) 



( 21 ) 



20 



Using (18) and (19) between Regions II and III gives the 
final equation in terms of m and a only: 



Because there is no way to solve (22) analytically, it 
is necessary to find the roots numerically using the IBM 
360/67 computer system at the Naval Postgraduate School. 

For a fixed iC = cr/f a search routine was used to find the 
several m's satisfying the equation. The computer work is 
described in Appendix A. 





( 22 ), 
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RESULTS AND CONCLUSIONS 



A . REVIEW 

There are a number of questions to be answered. MYSAK 
(1968a) in his finite-width model founds 

1. Shelf-wave numbers are inversely related to the 
shelf width, for a fixed frequency. 

2. There is a low wave number cut-off for edgewaves 
which is a function of the shelf width. That is, as the 
shelf width increases, the smallest possible wave number 
decreases (the largest possible wavelength increases). 

3. The fundamental mode edgewaves of REID (1958) 
with periods greater than one pendulum day do not exist 
over a continental shelf of finite width. 

It should be noted that Mysak's solution is only ap- 
proximate in that the maximum shelf depth h is assumed much 
less than the ocean depth, D, leading to an approximation 
of the equation expressing continuity at the edge of the 
shelf (Mysak's equation (10)). His results (labeled MA 
below) depend on this approximation. Exact solutions of 
Mysak's equation (6) (corresponding to equation (22) in this 
work) , were also generated so that comparisons could be 
made among MA, an exact solution for Mysak's model (ME) and 
results of the two-slope model studied in Section II (TS). 
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B. COMPARISON OF MYSAK'S APPROXIMATE SOLUTION WITH AN 
EXACT SOLUTION FOR MYSAK'S MODEL 

Nine cases were studied in order to compare ME and MA 
Two of the cases are illustrated in Fig. 2. In all cases 
the continental shelf is 100 km wide with a slope of .002,, 
duplicating Mysak's sample calculation. The bottom depth 
varied from 500m to 5000m. The shelf wave results are 
shown in Figs. 5-9. 



Note that: 

1. The approximate solution consistently gives a 
smaller wave number for a particular 05 and mode. It ap- 
pears to be the limiting condition for the exact solution. 

2. Except for the fundamental mode, an error of less 
than 1 % exists between corresponding modes of MA and ME in 
cases where the depth ratio h/D is smaller than .067. 

3. For the fundamental mode there . is still a signifi 
cant error introduced in m when using MA for large depth 
ratios and frequencies (>10% for 05 > 0.8 when D = 5 000m) . 




OCQO rn 






(Vertical exaggeration = 100:1) 



Fig. 2. Profile of the finite-width model 
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4. The edgewave results are shown in Figs. 18 - 20. 
Neither MA nor ME gives a fundamental edgewave mode similar 
to that of Reid. This can be seen directly from Mysak's 
equation (10) in the approximate case. Because the argu- 
ment is positive by definition, a must be negative in order 
for the Laguerre function to have roots (zeroes) . This in 
turn requires that either f > O’ or cr » f. 

5. A low wave number cut-off does exist for ME edge- 
waves, which diminishes with an increase in depth and in- 
creases with an increase of | CO | (Table I) . The cut-offs 
are always lower than those for MA, and are quite symmetric 
with respect to direction of travel. 



Table I. ME edgewave wave number cut-off (mW) 



Deep-water depth 
1000 
3000 
5000 



(m) First mode 
0.58 
0.32 
0.26 



Second mode 
1.32 
0.76 
0.58 



As pointed out by Mysak and Reid, the edgewave disper- 
sion relation is not symmetrical with respect to direction 
of travel, due to the influence of the Coriolis parameter. 
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C. COMPARISON OF THE MYSAK MODEL WITH THE TWO-SLOPE MODEL 





(Vertical exaggeration = 100:1) 

Fig. 3. Two-slope model. 

The Mysak model was compared with a two-slope model 
whose continental slope was .05 (Fig. 3). The results are 
shown in Figs. 5-9. Note that: 

1. Except for the fundamental mode for small a,' , there 
is little similarity between the dispersion relations for 
the two models, especially for the large deep-water depths 
(i.e., wide continental slopes). For a deep-water depth of 
5000 m, both mode 2 and 3 waves of TS have smaller wave 
numbers than mode 2 of ME for any given frequency. 

2. The fundamental TS shelf wave does not asympto- 
tically approach CO = 1.0 for large wave numbers as does the 
corresponding ME wave. In fact, the fundamental wave now 
behaves like the fundamental edgewave of Reid. 

3. A low wave number cutoff is still present for 
edgewaves, but at a significantly lower wave number (for 
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mode 1, nearly 


half that 


of ME) . 


The cutoffs 


are no 


longer 


symmetric with 


respec t to 


direction of travel 


(Table 


II) . 


Table II 


Two-slope 


edgewave 


wave number 


cutoff 


(mW^) 


Deep-water depth oo < 


: 0 


U) > i 


0 




m 


Mode 1 


Mode 2 


Mode 1 


Mode 2 




1000 


0.42 


1.32 


0.28 


1.28 




3000 


0.26 


0.76 


0.16 


0.76 




5000 


0.22 


0.58 


0.12 


0.58 





COMPARISON OF CONTINENTAL SLOPES 





Slope = 0.03 



Slope = 0.08 



(Vertical exaggeration = 100:1) 



Fig. 4. Two-slope models with different continental 
slopes . 



Three values were used for the value of the continen- 
tal slope: .03, .05, and .08 (Fig. 4). The results are 

shown in Figs. 10 - 18, for deep-water depths from 500 m 
to 5000 m. 



26 



1. Wave numbers for a particular mode of trapped 
waves over a fixed deep-water depth decrease with an in- 
crease of continental slope width. 

2. For the fundamental mode over a constant gradient 
continental slope, wave numbers increase with an increase 
in slope width (i.e., an increase in deep-water depth). 

All other modes decrease with an increase in width. 

3. A curve is not available for mode 3 for a conti- 
nental slope of 0.03 and deep-water depth of 5000 m. This 
is attributed to unknown problems of the computer routine. 
This problem does not occur elsewhere. 

4. Discontinuities appear in the dispersion relations 
for modes 2 and 3 with a continental slope 0.05 (in the 
deep-water depth range 2200-3400 m) and with a continental 
slope 0.08 (for depths greater than 2800 m) , suggesting the 
presence of complex values of m. A similar phenomenon is 
not observed for a slope 0.03. Equation (22) was investi- 
gated for complex values of m and real CO (Appendix B) and 
complex roots were found for a slope of 0.05 and depth 
2800 m (Fig. 22) . Since the surface height is the real 
part of C (x, y, t) complex values of m imply a spatial 
growth rate exp Cpm(my) } in the positive y direction. The 
roots m are complex conjugates, so that one wave grows and 
one decays at this rate. Then the most unstable wave 
(i.e., the one with the maximum growth rate) would be ex- 
pected to dominate the she If -wave spectrum. 
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5. No complex wave numbers were found for the edge- 
waves studied. Continental slope width is inversely pro- 
portional to wave number as in Mysak's results. 

E. RECOMMENDED FURTHER STUDY 

The next step would be to study further the effect of 
different continental shelf slopes using the TS model. 

More study is mandatory in D 4 above, both to? 

1. find its limits and 

2. determine if this is a mathematical curiosity or 
a physical reality. 

Future investigators should seek to avoid approxima- 
tions to their models. As this study has shown, the results 
for the exact equations can be significantly different than 
the approximations. 
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TWO -SLOPE SHELF ROUES, DEEP WATER- 5QQ HETERS 
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Fig. 6 
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TWO-SLOPE SHELF WflUES, DEEP WRTER= 1000 METERS 
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TWO-SLOPE SHELF WRUES, DEEP WRTER= 2000 METERS 
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TWO-SLOPE SHELF WRUES,DEEP WRTER= 3500 METERS 



Fig. 9 
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TWO-SLOPE SHELF LlflUES, DEEP WflTER= 5000 METERS 




o 

w 



o 



q 

V0 



o 

to 



o 

* 



o 

ro 



q 



o 



o 

a 






34 



EFFECT OF CONTINENTAL SLOPE, DEEP WATER- 500 METERS 
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EFFECT OF CONTINENTAL SLOPE, DEEP NATER= 1000 METERS 
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EFFECT OF CONT1NEMTRL SLOPE, DEEP WRTER= 2000 METERS 
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EFFECT OF CONTINENTAL SLOPE, DEEP WATER= Z500 METERS 



Fig. 14 
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EFFECT OF CONTINENTAL SLOPE, DEEP WATER= 2800 METERS 
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EFFECT OF CONTINENTAL SLOPE, DEEP WATER= 3000 METERS 



Fig. 16 
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EFFECT OF CONTINENTAL SLOPE, DEEP WATER= 3250HETERS 
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EFFECT OF CONTINENTAL SLOPE, DEEP WATER- 3SOO METERS 
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EFFECT OF CONTINENTAL SLOPE, DEEP WflTER= 5000 METERS 
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Fig. 20 
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Fig. 22 
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SHELF AND EDGE WAVES 
DEPTH =5000 METERS, SLOPE=.05 
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APPENDIX A 

NUMERICAL SOLUTION OF ONE-AND TWO-SLOPE MODELS 
// EXEC FOR TCLGP,DARM.eCRT='LI ST, MAP* , REG I GN. GC=95K ,T I ME. GO= 
/ / FOR T . $ YS I N PC * 

R EA L*8 L,M,H,COR ,SIG,SLO(2) ,W(2) ,ANS ,X ( 3 ) , 7 ( 3 ) , P ( 2 ) , 
GRAV, DELTA 

COMMON SLC, M .L,H, COR, SIG,W,X,Z, GRAV, A,KCCL 
DIMENSION AF (2 ) ,A( 3) ,ANSI (3030) , OM ( 60 ) , PT ( 6U , 3 ) , KLUF(8 
DATA PT,K(.UE /48C*0. 0,8*0/ 

READ (5, 1C 2) KOMAND 

IF KOMAND EQUALS C ( R LANK CARD IN DA T A OFCK) THE MYSAK 
SOLUTION AND APPROXIMATION ARE C CMPL) T ED. IF KOMAND EQUALS 
ANY INTEGER ( THROUGH 099) THE 2-SLCPE SOLUTION IS COMPUTED.. 

IF( KOMAND. NE.O) GOTO 4 
CALL M Y C A K 
GOTO 5Ct 
4 GRA V= 9 . 8D C 2 
DELTA=C. 2000-08 
TW0PI=nARCQC(-1.0C)* 2. DO 
READ( 5,52) W,COP,SLO 
WRITE! 6 , 56 ) W ,C OP , $ LO 
89 WW=W( 1 ) +W ( 2) 

KOUNT=C 
K TP L = 1 
I M= 8 

H=W( 1 )*SLO( l) +W(2) **10(2) 

CR IT^TWCPI /H+0.4CGD-C1 
NUM=CR I T/DELTA 
WR ITE ( 6 , 62 ) 

KEY=0 

3 DO 2 K = 1 ,60 

M=DFLOAT( KTRL) *PELTA 
I K= 1 

C S = FLOAT( K ) / 6 . E 0 1 

SIG=-C R*COR 

OM(K)=CS 

DO 5 J=KTRL,NUM 

LOOK = C 

00 1 1=1,2 

P( I )=( $IG*SIG-CPP*COR) / ( GR A V*SL C ( I ) )-CQR*M/S IG 
A E ( I ) = ( M- P ( I ) ) /( 2.D0*M) 

1 A( I ) =A E ( I ) 

A ( 3 ) = A ( 2) 

CALL FTNSIG(ANS,KEY) 

I F ( KOCL.NE.C ) GOTO ICC 
G0 T 0 41 

9 DO 8 M7=J,5CC 

M7 1 = M 7— 1 

8 ANSI( M Z)=AN C T( M 71) 

GO T 0 5 

41 TF(ANS.LF.C.DG) LOOK=1 
ANSI ( J )=ANS 

21 IF(L.L p .C.9°D-19)ANSI ( J)=0.D3 
JJ=J-1 

IF( ANSI ( J)*ANSI ( JJJ.LE.O.DC.AND.J.GT.l ) GOTO 30 
GOTO 5 

3C D T( K, IK ) = ( ANSI ( J) /(ANSI ( JJ)-AN? I ( J) ) *0E LT A + M ) *WW 
I c ( IK.E0.1.AND.J.GE.2) KTRL=J-1 
WPITE (6,61) P T ( K , I K ) 

I F ( I K .EQ . I v ) GOTO 31 
I K= I K + 1 
*5 M = M +DE L T A 

60 FORMA T ( 1H , 2HM = ,C1 C . 3 ,2HT0 , DIO . 3 ,40 X , l c HCOR I GL IS /S I GM A 
*=,D10.3) 

ICC DO 32 I M R = IK,IM 
32 KLUF ( IMp )=K-1 
I M= I K— 1 

31 WRI^E (6,60) DELTA, M,CS 

2 CONTINUE 

PT( 6C , 1 ) =7 . 5DOO 
DO 44 LEA K = 1 , 1 M 
A 4 KLUE ( LEAK ) =6C 
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00 90 KK= 1 , 8 
MJ = KLl'F( KK ) 

WR I TE ( 7 , 7 1 ) H.KK.MJ 

WRITE( 7,72) ( PTU I ,KK) ,11 = 1, «J) 

C C WR I TF ( 6,92) V J , ( PT ( I I , KK ) , I I =1 , M J) 

WR!TF( 6,62) 

CALL PLrTPl PT( 1,1) ,0M ,KLIJE ( 1 ) ,1 ) 

DO 33 IFK=2,8 

!F( KLUF( I SK) .LF.4) 0070 34 

33 CONTINUE 

34 I SK 1= ! SK- 1 

DO 35 kric=2,IFK1 

35 CALL PLCTF ( o j ( i , KRI C ) ,QM , K L'JE ( KRT C ) ,2 ) 

CALL PLOTF f PT« 1 , 1 SK ) , QM , KLUE ( I SK ) , 3 ) 

WRITE (6,65) 

R EA 0 ( 5 , P p ) W f 2 ) 

WR I TE ( 6 , 5 6 ) W ,C OR » FLO 

1 F ( W( 2 ) .NE .C.DO) GOTO 89 
5CC STOP 

52 FORMA T( « 1C . 3) 
c 3 FORMA T( 5E 1C. 3) 

64 FORMA T ( 1H , 5D2C.7) 

55 FOP MA T ( lH0,12X,lHM t l9X,lHL,17X,5HS!GMA t l7X,6HANFWER, 
215X,1HF,/,5D20. 7) 

56 FORMAT! 1H0,5C2C. 7) 

57 FORMAT ( 1H ,3025.12) 

*a cqoma T( 1HC,8X,4HP( 1) ,16X,4HP(2) ,16X ,4HA ( 1) , 16X,4HA( 2) , 
416X,^HA( 3) , / , 502 C. 9) 

61 FORMAT! 17H MIN OCCURS AT M=,E12.4) 

62 FORMA T( 1H 1 ) 

63 FORMAT ( 1 BH ZERO OCCURS AT *=, = 12. 4) 

65 FORMAT! 1HC,35X,2HMW) 

71 FORMA T ( E 20 . 7 , 21 1C) 

72 FOPMAT (5E16.7) 

88 FORMAT (010.3) 

c 2 FORMAT! 1HC, I 5,( /.5E20.7) ) 

1C2 F0PMATU2) 

END 



FURR OU T I NE M VSA K 

THIS PROGRAM SOLVES FOR MYSAK'S SOLUTION AND APPROXIMATION. 
FOR 2-SLOPE RETURN TO T HE MAIN PROGRAM. 

REAL*P 0, DD, COP ,W,GPA V,DEL ,OMEG A,LAMDA,DOELT A, ARGO, ANS 
»F p R I , SOT,ANSK,ANSWER,OUANT,M,SIG f DELT A 
DIMENSION QM( 60 ) ,ANSI (600) ,ANSJ(60C ) , PT ( 60 , 8 ) , KLUE ( P ) , 
PTT( 60 ,8) ,KLU! « ) 

DATA PT,PT7,KLUF ,KLU/ 96 0*0. 0,16*0/ 

READ (5,52) D,DD,COR, W 
WR I TE ( 6 , 5 6 ) D ,DD ,C OR , W 
WRITE! 6,6 2) 

GR A V= 9 . 8D2 
KOUNT=C 
P5 I M=4 
K TR L = 1 

3 DO 2 K=1 , 55 

CS=FLOAT( K ) /6.E1 

M=DFLOAT( KTRL) *0. 2 COD-O 8 

DFLTA=C.200D-08 

WR I T F ( 6,60 DELTA,M,CS 

SIG=CS*CCP 

I K= l 

I K K = 1 

DO 5 J=KTRL,500 
JJ=J 

IF! JJ.EO. 500. AND. I K • LF . I M ) GOTO 3C 
IF! JJ.F0.500.AND.IKK.LE. IM) GOTO 31 
L0OK=Q 
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dell=oelta 
Q M ( K ) - C S 
KEY = 0 

A=C.RD0-(C0R /SIG+( SIG*SIG-COR*CCR) / (GRAV*D/W4M) J/2.C0 

del=d /nr 

doelta= r.CR*cop*v<*w/(GPAv*n) 

lamda = m*w 

OMFGA = S IG /C OR 

ARGU=2.DC*LAVDA 

CALL OLA P ( AN S » A RGU ,A ,NN) 

C * LL ORDLA.R( A .ARGUtFPRI) 

OUANT=l.cr:+(DDELTA*DEL*(l .ODO- (OMFGA*OMFGA) ) / ( L AMDA* 
LAMOA ) 

SQT = nSQPT(QL’ANT) 

ANSWER =( $0 T - ( l.GDO /OMEGA ) -DEL* ( 1.0 DO -1.0 DO/ OMEGA ) )*ANS 
+ 2.CDJ*DEL*FPRI 
16 ANSI ( J )=ANS 
IF ANS J( J ) =A N SWEP 

IFdK.GT. IW.ANO.IKK.GT.IM) GOTO 2 
JJJ=J J-l 

IF( J.GF.2.AN0.ANSI ( JJ)*ANSI ( J J J ) . LF .0 . DO > GOT C 11 
I F( J .GF . 2 • AND. A NSJ ( J J ) *AN$ J( J J J ) . LE .0 . DO ) GGTC 12 
GOTO 5 

11 DI C =M+(ANST ( JJ) /(ANM ( JJJ) -ANSI ( JJ) )*OFLTA) 

IF( IK.FQ. l.ANO. JJ.GT.2) KTRL=JJ-2 

PT( K, IK )= OIF 
IK= I K + 1 

DIF1=M+(ANSJ( JJ) /(ANS J( JJJ)-ANE J( JJ) ) ♦DELTA) 

IF( ANSJ( JJ) ♦ANSJf J JJ) .GT.O.OD ) GOTO 13 
PTT(K f IKK)= DIF1 

ikk=ikk+i 

WRITE(6,65) PIFtDIFl 
GOTO 1 

13 WRITE(6,64) 0I C 
GOTO 1 

12 0IF1=M+(ANSJ(JJ)/(ANSJ( J J J ) - A NS J ( J J ) ) ♦DELTA) 
PTT(K»IKK)= 0 I F 1 

I KK = I K K + 1 

WR I TF ( 6 » 6 6 ) DIF1 

1 I F < IK.GT.Iw.tND.IKK.GT.IV) GOTO 2 

5 M=M +0 E L T A 
GOTO 2 

3C 00 32 I MF = I K , IM 

32 KLUF ( IMP )=K-1 

IF( JJ.FO.500.AMD.IKK.LE.IM) G0 T 0 31 
GOTO 3 A 

31 00 33 IMK=IKK, IV 

33 KLU( I M K ) =K— 1 

34 IF( IK.LT.IKK) IKK=IK 
! M= I KK- 1 

2 CONTI NL'F 

DO 14 K I N= 1 , 1 M 
KLU( K I N ) = 5 9 

14 KL'JFt KIN) = 59 

DO 105 M?7=l t 4 

MJ = KLL'F(M77) 

WRITE( 7,7CC)0D,MZZ ,MJ 

WR I TF ( 7 » 7G 1 ) (PT ( I I , MZZ ) , 1 1 =1 , VJ) 

M J = K L U ( MZZ ) 

WR I TE ( 7, 700 ) DO.MZZ.MJ 

WRITF{ 7,701) (PTT(II,MZZ) *11=1, MJ) 

DO 3G0 K°R = 1 ,6C 
PT( KPR,WZZ) = p T( KPR , MZ Z ) ♦ W 
3CQ PTT(KPR .MZZ) =PTT(KPR ,MZZ)^W 
M J = KL(JE ( M Z 7 ) 

WR I tc ( 7 , 7CC ) 0D,MZ7 t MJ 

WR I TE ( 7,701) (P*- (II , MZZ) ,11=1, WJ) 

M J = K L U (MZZ) 

WR I TE ( 7,700 0D,MZ7,MJ 

WRITE (7,701) ( P TT ( I I ,MZZ) ,11=1, MJ) 

105 CONTINUF 

6 DO °0 KK=1,S 
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MJ=KLUF( KK ) 

MJ.) = KLU(KK) 

WRITE! 6,92) MJJ. ( PTT ( I I ,KK) • I 1 = 1 , M J J ) 

9C WR I TF ( 6,92) MJ , ( PT ( I I , KK ) , 1 1 =1 , M J ) 

WRITF! 6,71 ) 

CALL PLOTP(PTT( 1 ,1 ) ,<JM,KLU(1 ) ,1 ) 

DO 35 MLK=2,7 

CALL PLCTP(PTT( 1 ,MLK) ,QM , K LU ( ML K ) ,2 ) 

35 CONTINUE 

CALL °l. CTP(PTT(1,8) ,QM,KLU(8) ,3) 

WRITF( 6,72) 

WR T TF ( 6,73) 

CALL PLCTP(PT(1 ,1) ,QM,KLUE(1) ,1 ) 

DO 36 Mpp=2,7 

CALL PLCTP! PT! 1 ,MOP) ,QM , KLIIE ( MOP ) , 2 ) 

36 CONTI NUF 

CALL PLCTP(PT(1,8) ,QP,KLUE<8),3) 

WR I TE ( 6,72) 

R EAD ( 5,65) DO 
WR I T c ( 6,56) D ,DD ,C OR , W 
IF! DD.NF.C.DC) OOTO 59 
RE T U D N 

c 2 F0PMAT(4F13. 3) 

5 3 FORMA T ( 5F 10. 3 ) 

54 FORMAT ( 1H ,5D2C.7) 

5 C FOP MAT ( 1HO,I2X,1HM,1PX,1HL,1 7 X,5HSIOMA,I7X,6HANSWFR, 15 
*5D2t . 7 ) 

56 FORMAT! 1HO ,4020. 7) 

57 FOPMATdH ,3025.12) 

58 FORMA T( 1HC,8X,4HP( 1) ,16X,4HP(2) , 16 X ,4HA ( 1 ) , 16X , AHA ( 2 ) , 
A /, 502C .9 ) 

60 F0pm*T( 3H M =»DlC.3»2HTC,DlJ.3»40X»15HSIGMA/CQRIQLI$=» 
DIG. 3) 

62 FORMAT(lHl) 

6 A FORMAT! 35H THE APPROXIMATION HAS A ROOT AT M=,F15.6,40 
H • T H5 EQUATION COES NOT WAVE A ROOT HERE.) 

65 FOR MA T ( 35H THE APPROXIMATION HAS A ROOT AT M=,E15.6,27 

H . THE EQUATION HAS A ROOT AT,S15.6) 

66 q 0P MA T ( 27H THE EQUATION HAS A RCCT AT , E 15 . 6 , 49H .THE AP 

PROXI MATI ON OOFS NOT HAVE A SOLUTION HERE.) 

71 FORMAT! 1H 1 ,2 7X » 2 6H PLOT OF MYSAK 1 $ ENTIRE EQN ) 

72 FQPMA T( 1WC,4CX,2HMW) 

73 FORMAT! 1H1,25X,29HPL0T OF MYSAK'S APPROX IMAT ICN ) 
ee FORMA T ( P 10. 3 ) 

92 FORMAT! 1HC, 15, ( /,5E2Q.d ) 

7CC FORMA T(F2C.7,2I 1C) 

7C1 FORMA T( 5F 12.4) 

END 



SUBROUTINE FINSIG!AN c ,KEY) 

REALMS L,M,X!3),F(3),FPRI(3),G(3),GPRI ( 3 ) , C , H , GR AV , COR 
» SIG, W( 2 ) » SLO! 2) ,ANF .7 (3) 

COMMON SL0,M,L,H,C0P,SIG,W,X,Z,GRAV, a,kocl 
DIMENSION A! 3) 

KOOL =0 

C = ( M*M*GP A V*H+C CR*COP-$IG*SIG)/ ( GR A V*H ) 

IF (C. . LF.O.OOGOTO 1 1 
L = DSQR T ( C ) 

X! 1 ) = W( 1 ) 

X! 2 ) = SLC( 1)*W(1) /SLO! 2 ) 

X( 3 ) = X ( 2 ) +W( 2) 

DO 2 1=1,3 

Z( I )=X( I )*2.DC*M 

CALL DLAP!F( I) ,7(1 ) ,A(I ) ,NN) 

CALL DLAPGG(F(I ) ,7(1 ) ,G(I) .A! I) ) 

CALL DR0LAP(A(I ) ,Z(I) ,FPRI (I) ) 

CALL DGDLA P ( A ( I ) , Z ( I ) , F op I ( I ) , G ( I ) , G PR I ( I ) ) 

IF(F( I ).GE.D,1D3C.0R.G(I).GF.0. 1030.0R.FPRI (I).GE. 

50. 103C.0R .GPRI ( I J.GE.G.1D30) GUTO ' 
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2 CONTINUE 

4N«=(L-«)*{F(3I*<F(I)*GPPI (2)-G(2)*FPRI Cl) )+G(3)*(F(2) 
2*FPR I ( 1 )-F( 1 )*FRRI (2) ) H-2*M*(FPni (3) *< FU)*GPRI (2 )- 
3G( 2)*FPR I( I n + GPPI ( 3) *(F(I) *FPRI ( 1) -F ( 1 ) *FPp I (2 ) ) ) 
TF(KEY.FO.l) WR I TE ( 6 , 5 1 ) ANS,M,L 
IF(AN$.GE.C. 1DOM GOTC 4 
I Rptijrm 
II L=9.99D-2C 
goto i 
4 KOOL = I 

WRITE (6,61) V 

r,o r n i 

51 Fnp MAT ( l*- 1 , 3D2C . 7 ) 

61 F0PM4 T( 34HGSGLUTI0N DCFS NCT EXIST BEYONT M=,£12.4) 

END 



SUBROUTINE 0 LA P f Y , X , A , NN) 

DlfF SOLVES FOP LAGIJFRPE FUNCTIONS OF 7 HE °IPST KING. GENER- 
ALLY thF ARGUMENT, A , WOULD RE EXPECTED TO BE NEGAT I V E. W HEN A 
IS PO$I T IVE ,NN IS SET TO 1 . 



REAL* 0 Y, X ,YY, YX,Y* ,YK ,YN 
YN=A 

YY= 1 , DC 
YK = 1 , QC 
Y X = I . D 0 
YZ= 1 ♦ DC 
I F ( A ) 3,1,2 

2 NN= 1 

3 YX=X/YZ*YN/YZ*YX 
YY=YY+ YX 

YN= YN+ 1 . DO 
YZ = Y7 + l.DC 

IF(DABS( YX) .GT.0.50E-C7) GOTO 3 
1 Y = YY 
RETURN 

EC FORMA T ( IH , 3 ( 5 X , E 1 4. 7 ) ) 

El FORMAT ( 1H , 5D2G. 7 ) 

END 



SUBROUTINE OLA PGG ( Y, X ,0 L , A ) 

DLAPGG SOLVES EG 0 LAGUEPRE FUNCTIONS OF THE SECOND KIND. WHEN 
A APPROACHES A NEGATIVF I MTEGEP , CCNTR CL IS TRANSFERRED TO 
DL AF G. WHEN CONVERGENCE DOE c NOT T AKE °LAC E BEFORE AN ANSWER 
IS OBTAINED, CONTROL IS RETURNED TO THE CALLING ROUTINE AND A 
MFSS4G C PRINTFDjWHEN X=C , THIS FUNCTION I c I NDFT C R MINATE. 

PEAL* 0 DN,D0,DX,D2,X,Y,DL,DA 
K0UNT=0 

I F ( X.FO.C ) GOTO 6 
DN= 1 .DC 
00=1. OC 
D X= l . DC 
DL=Y*DLOG( X) 

DA=A 
D2=0. DC 

4 DX=DX*DA /ON*X/DN 

D2=D2 + ( 1 . DC /DA ) - ( 2 . DO /DN) 

DQ=DX*D 2 
DL=DL +DP 
DA = DA + 1 . DC 

1/DA APPROACHING INFINITY? 
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IF!DABS!DA).LE.C. 10-02) GOTO 1 



ANSWER TOO LAPGF WITHOUT CONVERGENCE? 

IF <OABS(DX).GE.l.D65. OR. OABS!DX).LE. 1.0-65) GOTO IC 
CONVERGENCE? 

IF(DABS(DQ).LE.C.5D-C8. AND. KOUNT ,G T . 2 > GOTO 3 

DN=DN+l.DO 

KOUNT = KOUN T + I 

G0 T 0 4 

1 CALL DLA PG ( Y , X ,DL , A ) 

GOTO 3 

2 WRITE! 6,54)KCUNT 

3 RETURN 

10 WRITE<6,57) KOUNT, DX,DL 
GOTO 2 

6 WP I TF ( 6 * 5 3 ) 

GOTO 3 

53 FOR MA T ( 5 EH X IS C* LA GUERRE FUNCTION OF THE SECONC KIND 

DOFS NCT F XI ST) 

54 FORMAT! 1H ,1C0X,I3) 

K 5 FORMAT! 1H , 602C . ? , / , 02 C •. 7 ) 

E7 FORMA T( 2FH OVERFLOW A P PRO ACH I NG . AFT FR t I 3 , 1 5 H ITER- 
ATIONS ,D X= ,F 12 . 5 , RH ANO DL = » E12. 5 ) 

END 



SUBROUTINE DLAPG!Y,X,OL, A) 

OLAPG SOLVES FOR LAGUERPE FUNCTIONS OF THE SECOND KIND WHEN 
A IS A NEGATIVE IN T EGER. TO AVOID DIVIDING BY ZERO, THE NTH 
TFPM IS THF SUM OF N- 1 MULTIPLICATIONS. 

P EA L *8 Y, X,DA,77 ,DZ ,D1 ,02 ,0N,0L ,011 , DX ,D0 

KOUN T = 0 

DN= 1 .DC 

DQ= 1 . DC 

D X= 1 .DC 

DL=Y*DLOG( X) 

OA=A 
DZ= 1 . DC 
D2=C .DC 
22 DX=OX*X 
Dii=c.nc 
DO 1 I=C, KOUNT 
D 1= 1 . DC 
K0CK=I+1 
DO 2 K=C , I 
IF! K .EQ. I) GOTO 4 
IF(DABS!D1) .LE.C.50-1G) GOTO 1 

2 0 1=D 1 * ( DA +DF LOA T ( K ) ) 

A IF! KOOK.GT.KCUNT) GOTO 1 
pn 3 L=KOCK, KOUNT 
IF!DABS(D1).LE.C.5D-1C) GOTO 1 

3 D1=D1*(DA+DFL0AT(L) l 
1 011=0 11+01 

D2=D2-< 2.00/DN) 

OQ=DQ *DN*DN 

IF! DO. GE.C.1D40) GOTO 56 

ZZ=( D1 1 +D 1*! DA +DF LOA T( KOUNT) ) *021 *DX/0Q 

K0UNT=K0UNT+1 

I p( KOUNT .GT. 25) GOTO 56 

0L=DL + Z7 

IF(DABS!7Z ) .LE.C.5D-C8) GOTO 21 
0N = DN + 1 . DC 
GOTO 22 
21 RETURN 
56 WP I TF ! 6 , 53 ) 71 
G0 T 0 21 
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51 FORMAT! 1H ,M5X,E15.8) ) 

53 FOPMAT(36H G ( X) DID NCT CONVERT. LAST TFRM WAS,E15.9) 
END 



SUBROUTINE DRD LA P ( A , X , R F SULT ) 

DRDLAP SOLVES THE FIRST DERIVATIVE OF THE LAGUFRRF FUNCTION 
OF THE FIRST KIND. 

REAL* 0 YDA,YDD, RESULT, YDB,YDN,X 

yoa=a + i.dc 
rfsult=a 
YDN=1.D0 
YDP-6 

22 YDB = YDN+ l .DO 
YDD=YDD*YDA*X/ ( YDN*YDB ) 

RFSULT=PF<:ULT+YDD 

IP(DABS( YDD) .LT.0.5D-8) GOTO 21 
YDA= YDA + 1.00 
YDN=YDN+1 .DO 
GOTO 22 
21 RE T UPN 

23 FOR MA T ( 1H ,2 ( 5X ,E 1 5. Q ) > 

5 A FORMAT ( 1H ,6020.7) 

END 



SUBROUTINE DGD LA P ( A , X , RE SULT , PLG , RES t»L > 

DGOLAP SOLVES THE FID$T DERIVATIVE OF THE SECOND KIND CF 
LAGUERPE FUNCTION. WHEN X=u, THIS FUNCTION DOFS NOT EXIST. 

P EA L * 0 X, RES UL, RESULT, XA, PLG, XN,X3,X A, X2,XX t XZ,XC 
IFCX.EO.C) GOTO 1C 
XN= 2 . DO 
XA=A 

R E S UL =RE S ULT* 0 LOG ( X ) + P LG/ X + 1 . DO- ( 2 . DO *X A > 
X2=1.DC/X4-2.0C 
X X= XA 

21 XA = XA + 1 . DC 

IF(DABS(XA).LF.C.5D-''*) GCTO 22 

XX=XX*X/( XN*XN-XN) *XA 

X2= X2 -2.DC/XN +1.D0/XA 

XZ= XX* X2 

RFSUL=PFSUL + X7 

IF(DA8S< X7).LE.O.5D-0> GCTO 10 
XN= XN ♦ 1 . PC 

IF (DABS( X7).GE.1.D65) GCTO 23 
GO TQ 21 

22 KISS=1 

XN= 2 .DC 
XA =A 

RESUL=RESUL T *DLOG( X) + »LG/X +1 . DO - ( 2 . DC *X A ) 

X0=1 .DC 
XX= 1 . DC 
X2=- 2 . DC 
5 XX= X X* X 
X3=0. DC 

DO 6 11=0, KISS 
X A= 1 . DO 
K IN = K ISS+1 
DO 7 J.I=C,I! 

I F ( J J .EO. 1 1 ) GOTO 9 
IF(DA8S( XA).LE.G.5D-10) GOTO 6 
7 X4=XA*< XA+DFLOA T( JJ) ) 

9 IF ( KIN. GT. KISS) GOTO 6 
DO 8 LL=KIN,KISS 
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I F ( B S ( X4I.LE. 0.50-10) GOTO 6 
8 XA=X4*( XA+OFLOA T( LL) ) 

6 X 3= X3 + X4 

X2=X2-( 2.00/XN) 

XO=XO*XM*DFLCAT( KI SS) 

XZ=( X3+X4*( XA+OFLCATI KT SS) ) 4X2) *XX/XQ 

KISS=KISS+1 

RF$UL=PF$UU-X7 

IFC DAB S( X7 ) . LE . 0. 5D-8 ) GOTO 10 
XN=XN+ 1.00 

IF(KI$S.LT.25) GOTO c 
WR I TF ( 6,54) XZ 

23 N=XN 

WRITE! 6,24) N 
10 RFTUPN 

24 FORMA T ( **CH OVERFLOW ABOUT TO Of CUR IN DGDLAP AFT C R, 13, 
2 7H TERMC.) 

5 1 FORMAT! 1H , 4 ( 5 X , El 5 . 8 ) ) 

54 FORMA T ( 1H ,36HGMX) 010 NOT C CNVFRGF. LAST TERM WAS, 
1F15.8) 

END 

/ /GO . FT 06FGG 1 00 DCB = < PECF m = f A ,B LKS I 7 E=1 33 ) ,5 P ACE= ( CVL , ( 15 , 1 
//GO.SYSTN DO * 
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APPENDIX B 



NUMERICAL SOLUTION FOR TWO-SLOPE COMPLEX ROOTS 

// EXEC FORTCLGP, PARM. FORT= • L I ST, MAP* , REGION. G0=100K, TIME. GO 
//FORT .SYS IN DD * 



COMPLEX* 16 M,L,ANS,X( 3) ,Z(3) ,P(2) ,A(3) ,AE(2> .ANSI (1620 
REAL *8 H . COR . S IG, SLO( 2) , W ( 2 ) , GR AV , P AR 1 , PA R 2, ANSWER ( 162 
COMMON M,L»X»Z»A»H»COR»SIG»SLO»W»GRAV»KOOL 
DATA ANSI/1620*( O.EO.O.EO)/ 

THIS PROGRAM SOLVES FOR THE COMPLEX ROOTS OF THE TWO 
SLOPE MODEL. IN THIS CASE, THE ROOT WAS FOOND FOR OM- 
E GA= 12/60 AT J= 10 AND K= 30. THIS GIVES AN ANSWER OF 
( 1.5459r-07,0.1249E-07) FOR THE ROOT. 

GRA V=9. 8D2 

RE AD ( 5 , 52 ) W.COR.SLO 
WRITE(6,56)W»C0R» SLO 
99 WW=W ( 1 )+W( 2) 

H=W ( l ) *SLO ( 1 ) +W ( 2 ) *S LO ( 2 ) 

WRITE(6,100) 

PARl=1.541234E-07 
S IG=-12.00*COR/6.00E01 
DC 2 J=1 , 20 
PAR 1 = P AR 1 + 0 . 500D- 10 
M = P AR 1 * ( 1 . DO . 0 . DO ) 

DC 1 1=1,2 

P ( I )=( SI G*S I G-COR *COR ) / { GRAV* SLO ( I ) )-COR*M/SIG 
AE ( I ) = ( M- P ( I ) ) / ( 2 . D0*M ) 

1 All )=AE( I) 

A ( 3 ) =A ( 2 ) 

CALL F I NS I G ( ANS , KEY ) 

ANSI ( J)=ANS 

ANSWER (J )=CDABS( ANS I( J ) ) 

DC 2 K=22 1,260 

P AR2= DFLOAT ( K ) * C.50CD-10 

M=PAR1*( 1 .CO, 0.00 ) + P AR2 * (0.00,1. DO) 

DO 8 1=1,2 

P ( I )=( SIG*SIG-COR*COR)/(GRAV*SLO< I ) )-COR*M/SIG 
AE( I ) = ( M— P ( I ) ) / ( 2 . D0*M ) 

8 A ( I ) = AE( I ) 

A (3 ) = A ( 2 ) 

CALL FI NS IG( ANS , KEY ) 

KOFE=(K-22C>*20 + J 

PRINTS ANS I( 2 1 I - ANS I (820) 

ANSI ( KOF E ) = ANS 

ANSWER (KOF E) = CD AB S ( AN S I ( KOF E ) ) 

M = M— ( 2.DC*PAR2*(0.D0» l.DO) ) 

DO 9 1=1,2 

P( I)=( SIG*SIG-CCR*COR)/(GRAV*SLG< I) )-COR*M/SIG 
AE ( I )=(M-P( I ) )/( 2.D0*M) 

9 A( I ) = AE( I ) 

A( 3) =A( 2) 

CALL FINS IG( ANS, KEY) 

K0FF=(K-180) *20 + J 

PRINTS ANS I ( 821)- ANSK 1620) 

ANSI ( KOFF )=ANS 

ANSWER (KOFF)=CDABS( ANSI (KOFF) ) 

2 CONTINUE 

DC 4 LPA= 1,31 
K 1= ( LPA-1 ) *20+1 
K2=LPA*20 

WRITE (6, 101) ( ANSI (LP) ,ANSWER(LP) ,LP=K1,K2) 

4 CONTINUE 

52 FORMAT ( 5E 1 0. 3 ) 

56 FORMAT ( 1H0, 5D2C. 7 ) 

100 FORMAT(lHl) 

101 FORMAT( 1H0,( 9D14. 5) ) 

STOP 
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END 



SUBROUTINE F I NS I G < ANS , K EY ) 

COMPLEX* 16 M,L,X(3),Z(3),A(3) ,F < 3 ) ,FPR I < 3 ) ,G< 3) ,GPRI 
* ( 3 ) , C , ANS 

RF AL*8 H,COR,SIG,SLO( 2) ,W(2) ,GRAV, CRAPS 
COMMON M,L ,X,Z , A , H ,COR , S I G , SLO , W , GRA V ,KOOL 
K0CL=0 

CRAPS=0. ID 30 

C = ( M*M*GR A V*H+COR *COR-SIG*SIG)/(GRAV*H) 

L=CDSCRT(C) 

x< i)=wm 

X <2)=SLG< 1 )*W< 1 )/SLO( 2) 

X ( 3 ) =X ( 2 ) +to ( 2 ) 

DO 2 1=1.3 

Z ( I ) =X ( I ) *2 . DO *M 

CALL DLAP(F( I ) ,Z( I) ,A( I) ,NN) 

CALL DLAPGGIF < I ) ,Z( I ) ,G(I ),A( 1)1 
CALL DRDLAPI A ( I ) , Z(I) , FPR I ( I)) 

CALL DGDLAP ( A ( I ) , Z ( I ) ,FPRI < I ) , G ( I ) ,GPRI < I ) ) 

I F (CDARS ( F ( I ) ) .GE. CRAPS. OR. CDABSIGII )). GE .CRAPS. OR CD- 
ABSI FPRI ( I ) ) .GE. CRAPS. OR . CDAB S ( GPR I ( I ) ) .GE. CRAPS) GOTO 
? J 

2 CONTINUE 

ANS=(L-M)*(F(3 )*( F ( 1 ) *GPR I ( 2 ) -G ( 2 ) *FPR I ( 1 ) )*G ( 3 ) * ( F ( 2 ) 
2*F PR I ( 1 ) - F ( 1 ) *FRR 1(2) ))+2*M*(FPRI (3)*(F(1)*GPRI(2)- 
3G ( 2 ) *FPP I ( 1 ) ) +GPR I ( 3 ) * ( F ( I)*FPRI( 1 )-F ( 1 ) *FPR 1(2) ) ) 
IF(KEY.EO.l) WRITE(6,51) ANS,M,L 
IFICDABS(ANS) . GE.0.1D4) GOTO 4 
1 RETURN 
4 KOOL= 1 

WR I TE ( 6 » 6 1 ) M 
GOTO 1 

51 FORMAT ( 1H ,3D20.7) 

61 FORMAT I34H0S0LUT ICN DOES NOT EXIST BEYOND M=,E12.4> 

END 



SUBROUTINE DL AP ( Y , X , A , NN ) 

CLAP SOLVES FOR LAGUEPRE FUNCTIONS OF THE FIRST KIND. GENER- 
ALLY THE ARGUMENT, A, WOULD BE EXPECTED TO BE NEGAT I VE . WHEN A 
IS POSITIVE, NN IS SET TO 1. 

COMPLE X* 1 6 Y ,X,A,YN,YX 

RE AL *8 YY , YZ 

Y N=A 

Y Y= 1 . DO 

YK=1 .DO 

YX=1 . DO 

YZ=l.CO 

3 YX=X/YZ*YN/YZ*YX 
YY=YY+YX 
YN=YN+ 1 . DO 
YZ=YZ*1 • DO 

IFICDABS(YX) .GT.0.5D-07) GOTO 3 
1 Y= YY 
RETURN 

50 FORMAT ( 1H , 3 ( 5X , E 1 4 . 7 ) ) 

51 FORMAT (1H ,5D20.7) 

END 



SUBROUTINE DL APGG ( Y , X , DL , A ) 
CCMPLEX*16 Y, X,DL,A,DA,DQ,DX,D2 
REAL*8 ON 



57 



CLAPGG SOLVES FOR LAGUERRE FUNCTIONS OF THE SECOND KIND. WHEN 
A APPROACHES A NEGATIVE INT EGER , CONTROL IS TRANSFERRED TO 
OL APG. WHEN CONVERGENCE DOES NOT TAKE PLACE BEFORE AN ANSWER 
IS OBTAINED, CONTROL IS RETURNED TO THE CALLING ROUTINE AND A 
MESSAGE PRINTED;WHEN X=0, THIS FUNCTION IS INDETERMINATE. 

K0UNT=0 

IF (CDABS (X ) . EQ.O. CO) GOTO 6 

DN= 1 • DO 

DQ=1 .DO 

DX= I . DO 

DL = Y*C DLOG ( X ) 

D A=A 
D2=0 . DO 

A DX=DX*DA/DN*X/DN 

D2 = D2+ (1 .DO/DA )-( 2 .DO/DN) 

DO=DX*D2 
DL=DL+DU 
DA= D A+ 1 . DO 

1/DA APPROACHING INFINITY? 

IF (CDABSI CX) .GE.O . 1D65 .OR .CDABS ( DX ) .LE.O. ID-60) GOTO 1 
ANSWER TOO LARGE WITHOUT CONVERGENCE? 

IF (CDABS(DA) .LE.O. ID-2) GOTO 1 
CONVERGENCE? 

IFICDABSIDO) .LE.O. 5D- 8.AND.KOLNT.GT.2) GOTC 3 
DN=DN+ 1 . DO 
KOUNT=KOUNT+1 
GOTO A 

1 CALL CLAPGI Y,X,OL,A) 

GOTO 3 

2 WRITE! 6,5A)KOUNT 

3 RETURN 

10 WPITE(6»57) KCUN T , DX , CL 
GOTO 2 

6 W R I T E ( 6 , 53 ) 

GOTO 3 

53 FORMAT ( 59H X IS 0. LAGUERRE FUNCTION OF THE SECOND KIND 
a COES NOT FXIST) 

5A FORMAT ( IH ,100X,I3) 

55 FORMAT ( 1H , 6D20 . 7 , / , D 20 . 7 ) 

57 FORMAT ( 29H OVERFLOW APPROACH I NG . AFT ER ,I3,15H ITERAT- 
/ I0NS,DX=,E12. 5,9H AND DL=,E12.5) 

END 



SUBROUTINE OL APG ( Y , X, CL , A ) 

COMPLEX* 16 Y, X,A,DL»DA,ZZ ,D1»D11 ,DX 
RFAL*8 D2,DN,0Q 

DL APG SOLVES FOR LAGUERRE FUNCTIONS OF THE SECONC KIND WHEN 
A IS A NEGATIVE INTEGER. TO AVOID DIVIDING BY ZERO, THE NTH 
TERM IS THE SUM OF N-l MULTIPLICATIONS. 

KGUNT=0 
DN= 1 . DO 
DQ= 1 .DO 
DX=1 .DO 
DL=Y*C DLOG ( X ) 

DA=A 
DZ=1 .CO 
02=0. DO 
22 DX=DX*X 
Cl 1=0 .DO 
DO 1 I =0 , KOUNT 
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01=1. DO 
KCOK= I + 1 
002 K=0 » I 
IF(K .EQ. I ) GOTO A 
IFICDABS(Dl) .LE. 0.50-10) GOTO 1 

2 D1=D1*(DA+DFL0AT(K) ) 

A I F ( KOOK.GT .KUUNT ) GOTO 1 
DO 3 L=KOCK. KOUNT 
IFICDARSIDl) .LE. 0.50-10) GOTO 1 

3 01 = D1*(DA + DFLOAT(L) ) 

1 011=011+01 

02=02-1 2. DO/DN) 

DQ=DQ*DN*DN 

I F ( OQ . GE.0.10A0) GOTO 56 

Z Z= ( 0 1 1 + 01 *( DA +0 FLOAT (KOUNT) )*D2)*DX/DQ 
KCUNT=K0UNT+1 
I F ( KOUNT . G T . 2 5 ) GOTO 56 
DL =DL + ZZ 

IF(CDABS(ZZ) .LE.O .50-8) GOTO 21 
DN=DN+ 1.00 
GOTO 22 
21 RETURN 
56 WR I TE ( 6 » 5 3 ) ZZ 
GOTO 21 

51 FORMAT ( 1H , A ( 5X , E 1 5 . 8 ) I 

53 FORMAT ( 1H ,35HG(X) OID NOT CONVERGE . LAST TERM WAS,E15. 
ENO 



SUBROUTINE OROLAP l A, X , RESULT ) 

COMPL E X* 1 6 A, XfRESULT , YOA » YDD 
REAL*8 YOB, YON 

OROLAP SOLVES THE FIRST DERIVATIVE OF THE LAGUERRE FUNCTION 
OF THE FIRST KIND. 

Yf)A= A+ 1 . DC 
RESULT=A 
YCN=1 .00 
YDD= A 

22 YDB=YDN + 1 .00 
YCD=YDD*YDA*X/ (YCN*YDB) 

R ESULT=R E SUL T + YDO 

I F(CDABS(YDO) .LT. 0.50-8) GOTO 21 
YDA=YDA+1 .00 
Y CN= YDN+1 .00 
GOTO 22 
21 RETURN 

23 FORMAT ( 1H , 2 ( 5X , El 5 .8 ) ) 

5A FGRMAT ( 1H ,6020.7) 

END 



SUBROUTINE OGOLAP ( A, X , RESULT , PLG, RESUL ) 

COMPLEX* 16 A, X, RESULT , PLG , RESUL , X A , X 2 , XX , X Z , XA, X3 
RE AL*8 XN 

OGOLAP SOLVES THE FIRST DERIVATIVE OF THE SECOND KIND OF 
LAGUERRE FUNCT I ON . WHEN X=0, THIS FUNCTION OOES NOT EXIST. 

IF (CDABSIX ) .EQ.O.CO) GOTO 10 
XN=2 .00 
X A= A 

RESUL=RESULT *CDLOG (X) +PLG/X + 1 .00- ( 2 . DO*X A ) 
X2=l.DO/XA-2.DO 
XX = X A 

21 XA=XA+ 1 . DO 

IF(CDABSIXA) .LE.0.5D-A) GOTO 22 



59 



XX=XX*X/(XN*XN-XN)#XA 
X 2=X 2 -2.DC/XN +1.D0/XA 
XZ=XX*X2 
RESUL=RCSUL +X Z 

IFICDABSI XZ) .LE. 0.50-8) GOTO 10 
X N=X N + 1.D0 

IFICDABS(XZ) .GF.0.1D60) GOTO 23 
GO TO 21 

22 K I SS = 1 

XN=2 . DO 
X A=A 

RESUL=RFSULT*CDLOG(X) +PLG/X+1 .DO- I 2.D0*XA) 

XQ=1 .DO 
XX=1 .00 
X 2=-2 * DO 

5 XX=XX*X 
X 3=0 .00 

DO 6 I 1=0, KISS 
XA = 1 . DO 
K IN=KI SS+1 
DC 7 JJ=0, I I 
IF(JJ.EQ.II) GOTO 9 
IFICDABS(XA) .LE.G.5D-10) GOTO 6 

7 XA=XA*(XA+DFLOAT (JJ) ) 

9 IFIKIN.GT.KISS) GOTO 6 
OC 8 L L=K I N, K l SS 
IFICDABS(XA) .LE.0.5D-10) GOTO 6 

8 XA=X4*( XA+DFLOATI LL) ) 

6 X 3=X 3+XA 
X2=X2-(2.DC/XN) 

XC=XQ«XN*DFLOAT(KISS) 

XZ=(X3+XA*(XA+DFL0AT( KISS ) )*X2)*XX/XQ 

K I SS=K I SS+ 1 

RESUL=RESUL+XZ 

IFICDABS(XZ) .LE.0.5D-8) GOTO 10 
XN=XN+ 1 . DO 

IFIKISS.LT. 25) GOTO 5 
WRIT E ( 6, 5A ) XZ 

23 N=XN 

WP I TE ( 6 , 2 A ) N 
10 RETURN 

2 A FORMAT ( AOFI OVERFLOW ABOUT TO OCCUR IN DGDLAP AFTER, I 3 
+ 7H TERMS. ) 

51 FORMAT ( 1H ,A( 5X , E 15. 8 ) ) . 

5 A FORMAT ( 1 H ,36HG'(X) DID NOT CONVERGE. LAST TERM WAS, 

, E15.8) 

END 

// GC . FT 06F00 1 DC DC E= I R ECFM = FA , BLK S I ZE= 1 33 ) , SP ACE= ( CYL , ( 15, 1 
//GO. SYSUOUMP OD SYSOUT=A 
//GO.SYSIN DD * 

0.1 ODD 08 0.520D 07 0.729D-0A 0.2000-02 0.500D-01 
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